Projet télédétection / Master 2 ECOMONT

Ce projet est réalisé dans le cadre du Master 2 ECOMONT et porte sur l’étude spatiale des données SMOD (date de fonte de la neige) et NDVI (indice de végétation par différence normalisée) obtenues sur une zone délimitée, située sur le territoire de la communauté de communes de la vallée de Chamonix.

1 Chargement des librairies

library(raster)
library(dplyr)
library(mblm)
library(tiff)
library(rgdal)
library(rgeos)
library(terra)
library(basemaps)
library(mapedit)
library(rasterVis)
library(RColorBrewer)
library(tidyverse)

2 Sources des données

2.1 Données NDVI

Institution(s):
Pôle THEIA du CNES (traitement des images Sentinel-2)
CREA Mont-Blanc (calcul du NDVI)

Auteurs
Pôle THEIA - https://theia.cnes.fr/atdistrib/rocket/#/home
Brad Carlson (CREA Mont-Blanc)

2.2 Données CCVCMB

Institution(s): Communauté de Communes de la Vallée de Chamonix Mont-Blanc (CCVCMB)

2.3 Données SMOD

Données initiales
CESBIO Toulouse
Pôle THEIA du CNES
Accès : https://theia.cnes.fr/atdistrib/rocket/#/home
Citation : https://www.mdpi.com/2072-4292/12/18/2904/pdf
Contact :
Traitement des données :
CREA Mont Blanc

3 Objectifs du projet

Notre projet a consisté à étudier s’il existait une corrélation entre l’évolution du SMOD et l’évolution du verdissement observé (via le calcul du paramètre NDVI), entre 2017 et 2021, sur une zone délimitée de la communauté de communes de la vallée de Chamonix (CCVCMB).

4 Analyse spatiale et statistique des données

4.1 Import et description des couches raster

4.1.1 Import des couches NDVI par année

rastlist.2017.NDVI <- list.files(path = "S2_NDVI_CCVCMB", pattern='2017', all.files=T, full.names=T)
rasters.2017.NDVI <- stack(lapply(rastlist.2017.NDVI, raster))

rastlist.2018.NDVI <- list.files(path = "S2_NDVI_CCVCMB", pattern='2018', all.files=T, full.names=T)
rasters.2018.NDVI <- stack(lapply(rastlist.2018.NDVI, raster))

rastlist.2019.NDVI <- list.files(path = "S2_NDVI_CCVCMB", pattern='2019', all.files=TRUE, full.names=T)
rasters.2019.NDVI <- stack(lapply(rastlist.2019.NDVI, raster))

rastlist.2020.NDVI <- list.files(path = "S2_NDVI_CCVCMB", pattern='2020', all.files=TRUE, full.names=T)
rasters.2020.NDVI <- stack(lapply(rastlist.2020.NDVI, raster))

rastlist.2021.NDVI <- list.files(path = "S2_NDVI_CCVCMB", pattern='2021', all.files=TRUE, full.names=T)
rasters.2021.NDVI <- stack(lapply(rastlist.2021.NDVI, raster))

4.1.2 Caractéristiques des couches NDVI

## class      : RasterStack 
## dimensions : 2636, 2348, 6189328, 30  (nrow, ncol, ncell, nlayers)
## resolution : 10, 10  (x, y)
## extent     : 324930, 348410, 5077160, 5103520  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=32 +datum=WGS84 +units=m +no_defs 
## names      : SENTINEL2A_20200224_NDVI, SENTINEL2A_20200315_NDVI, SENTINEL2A_20200325_NDVI, SENTINEL2A_20200404_NDVI, SENTINEL2A_20200409_NDVI, SENTINEL2A_20200424_NDVI, SENTINEL2A_20200504_NDVI, SENTINEL2A_20200509_NDVI, SENTINEL2A_20200519_NDVI, SENTINEL2A_20200524_NDVI, SENTINEL2A_20200529_NDVI, SENTINEL2A_20200603_NDVI, SENTINEL2A_20200623_NDVI, SENTINEL2A_20200708_NDVI, SENTINEL2A_20200713_NDVI, ... 
## min values :               -0.9983560,               -0.9981877,               -0.9985891,               -0.9978938,               -0.9985730,               -0.9980031,               -0.9985617,               -0.9979624,               -0.9981035,               -0.9982855,               -0.9980518,               -0.9969445,               -0.9991574,               -0.9974279,               -0.9978781, ... 
## max values :                0.7155070,                0.7039071,                0.7296859,                0.8276377,                0.8230168,                0.8353869,                0.8483894,                0.9216663,                0.8404192,                0.8615274,                0.8427536,                0.8281259,                0.8500417,                0.8184949,                0.8514056, ...

Les rasters des données NDVI issues du satellite Sentinel2 ont une résolution de 10 m par 10 m et sont projetés dans le référentiel SCR WGS84 / UTM zone 32N.

4.1.3 Visualisation d’une couche NDVI (24/02/2020)

Figure 1: Couche NDVI limitée à la CCVCMB le 24/02/2020

Figure 1: Couche NDVI limitée à la CCVCMB le 24/02/2020

Les zones blanches correspondent à des nuages. Les couleurs les plus vertes dessinent la vallée de Chamonix et le gradient des couleurs jaunes à roses refléte le gradient altitudinal des versants qui encadrent la vallée.

4.2 Import des couches SMOD

4.2.1 Import des couches SMOD par année

SMOD_17 <- raster("Sentinel2_neige/SMOD_CCVCMB_2017.tif")
SMOD_18 <- raster("Sentinel2_neige/SMOD_CCVCMB_2018.tif")
SMOD_19 <- raster("Sentinel2_neige/SMOD_CCVCMB_2019.tif")
SMOD_20 <- raster("Sentinel2_neige/SMOD_CCVCMB_2020.tif")
SMOD_21 <- raster("Sentinel2_neige/SMOD_CCVCMB_2021.tif")

4.2.2 Caractéristiques des couches SMOD

## class      : RasterLayer 
## dimensions : 1318, 1175, 1548650  (nrow, ncol, ncell)
## resolution : 20, 20  (x, y)
## extent     : 324920, 348420, 5077160, 5103520  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=32 +datum=WGS84 +units=m +no_defs 
## source     : SMOD_CCVCMB_2017.tif 
## names      : SMOD_CCVCMB_2017 
## values     : 0, 243  (min, max)

Les rasters des données SMOD issues du satellite Sentinel2 ont une résolution de 20 m par 20 m et sont projetés dans le référentiel SCR WGS84 / UTM zone 32N. Il faut donc redimensionner les couches SMOD pour avoir la même résolution que les couches NDVI (10 m x 10 m).

4.2.3 Changement de la résolution des couches SMOD

SMOD_17 <- resample(x = SMOD_17, y = rasters.2020.NDVI[[1]], method = "bilinear")
SMOD_18 <- resample(x = SMOD_18, y = rasters.2020.NDVI[[1]], method = "bilinear")
SMOD_19 <- resample(x = SMOD_19, y = rasters.2020.NDVI[[1]], method = "bilinear")
SMOD_20 <- resample(x = SMOD_20, y = rasters.2020.NDVI[[1]], method = "bilinear")
SMOD_21 <- resample(x = SMOD_21, y = rasters.2020.NDVI[[1]], method = "bilinear")
SMOD_17
## class      : RasterLayer 
## dimensions : 2636, 2348, 6189328  (nrow, ncol, ncell)
## resolution : 10, 10  (x, y)
## extent     : 324930, 348410, 5077160, 5103520  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=32 +datum=WGS84 +units=m +no_defs 
## source     : memory
## names      : SMOD_CCVCMB_2017 
## values     : -22.25, 287  (min, max)

La résolution des couches SMOD a bien été modifiée.

4.2.4 Visualisation d’une couche SMOD (2017)

Figure 2: Couche SMOD en 2017 autour de la vallée de Chamonix

Figure 2: Couche SMOD en 2017 autour de la vallée de Chamonix

Des couleurs claires roses/jaunes sont associées à la vallée de Chamonix et donc à des zones pour lesquelles la neige fond précocement. Les couleurs vertes vives sont associées à des zones hautes en altitude pour lesquelles la neige fond tardivement ou jamais.

4.2.5 Application du masque CCVCMB aux couches SMOD

4.2.5.1 Import des limites de la communautés de communes (CCVCMB)

La couche correspondant aux limites de la CCVCMB est importée au préalable pour appliquer son emprise sur les couches SMOD afin d’obtenir des limites similaires aux couches NDVI.

CCVCMB <- shapefile("Communes/CC_Vallee_CHAMONIX.shp")
CCVCMB <- spTransform(CCVCMB, CRS = crs("+proj=utm +zone=32 +datum=WGS84 +units=m +no_defs "))

La couche CCVCMB est également reprojetée dans le même référentiel SCR que les couches NDVI.

4.2.5.2 Application du masque CCVCMB

SMOD_17M <- mask(SMOD_17, CCVCMB)
SMOD_18M<- mask(SMOD_18, CCVCMB)
SMOD_19M<- mask(SMOD_19, CCVCMB)
SMOD_20M<- mask(SMOD_20, CCVCMB)
SMOD_21M<- mask(SMOD_21, CCVCMB)

4.2.5.3 Visualisation de l’effet du masque CCVCMB

Figure 3: Visualisation de la couche SMOD 2017 après application du masque CCVCMBFigure 3: Visualisation de la couche SMOD 2017 après application du masque CCVCMB

Figure 3: Visualisation de la couche SMOD 2017 après application du masque CCVCMB

4.3 Préparation des tableaux NDVI max et SMOD

L’objectif est désormais de calculer le NDVI max pour chaque année afin de pouvoir calculer la pente liée à l’évolution du NDVI max sur les cinq années d’étude (2017-2021). La pente associée à l’évolution des couches SMOD sur les cinq mêmes années sera également calculée.

4.3.1 Calcul NDVI max par année

Les rasters et les tableaux associés aux NDVI max sont générés.

NDVImax.2017 <- max(rasters.2017.NDVI, na.rm=T)
NDVImax.2017_df <- as.matrix(NDVImax.2017)

NDVImax.2018 <- max(rasters.2018.NDVI, na.rm=T)
NDVImax.2018_df <- as.matrix(NDVImax.2018)

NDVImax.2019 <- max(rasters.2019.NDVI, na.rm=T)
NDVImax.2019_df <- as.matrix(NDVImax.2019)

NDVImax.2020 <- max(rasters.2020.NDVI, na.rm=T)
NDVImax.2020_df <- as.matrix(NDVImax.2020)

NDVImax.2021 <- max(rasters.2021.NDVI, na.rm=T)
NDVImax.2021_df <- as.matrix(NDVImax.2021)

4.3.2 Visualisation NDVI max par année

Figure 4: Couches NDVI max par annéeFigure 4: Couches NDVI max par année

Figure 4: Couches NDVI max par année

Peu de variations sont observées sur les 5 années. Un léger front de verdissement est observé sur le versant est de Chamonix mais il reste compliqué d’identifier des changements facilement observables.

4.3.3 Homogénéisation du nombre de pixels entre les 5 années pour les couches SMOD et NDVI max

Les masques des couches présentant les nombres de pixels différents de “NA” les plus faibles sont appliqués à l’ensemble des couches afin d’obtenir un nombre de pixels similaire pour l’ensemble des couches.

NDVImax.2017<-mask(NDVImax.2017, NDVImax.2020)
NDVImax.2018<-mask(NDVImax.2018, NDVImax.2020)
NDVImax.2019<-mask(NDVImax.2019, NDVImax.2020)
NDVImax.2020<-mask(NDVImax.2020, NDVImax.2020)
NDVImax.2021<-mask(NDVImax.2021, NDVImax.2020)
SMOD_17M <- mask(SMOD_17M, SMOD_18M)
SMOD_18M<- mask(SMOD_18M, SMOD_18M)
SMOD_19M<- mask(SMOD_19M, SMOD_18M)
SMOD_20M<- mask(SMOD_20M, SMOD_18M)
SMOD_21M<- mask(SMOD_21M, SMOD_18M)

NDVImax.2017 <- mask(NDVImax.2017, SMOD_17M)
NDVImax.2018 <- mask(NDVImax.2018, SMOD_17M)
NDVImax.2019 <- mask(NDVImax.2019, SMOD_17M)
NDVImax.2020 <- mask(NDVImax.2020, SMOD_17M)
NDVImax.2021 <- mask(NDVImax.2021, SMOD_17M)

Nombre de pixels des couches SMOD:

## [1] 3501381 3501381 3501381 3501381 3501381

Nombre de pixels des couches NDVI:

## [1] 3501381 3501381 3501381 3501381 3501381

Le nombre de pixels est désormais identique pour l’ensemble des couches NDVI max et SMOD.

4.3.4 Réduction de la zone d’étude

Le nombre de pixels étant très important, le script requiert 6 heures de calcul pour être exécuté. Dans le cadre de notre projet étudiant, le choix a donc été de délimiter une zone d’étude restreinte afin de réduire le temps de calcul nécessaire.

4.3.4.1 Extraction d’une portion des jeux de données raster

Une zone a été définie manuellement pour couvrir à la fois la vallée et les versants et a été appliquée aux rasters NDVI et SMOD, afin d’extraire une portion de ces rasters selon un modèle d’étendue.

cropbox2 <-c(336000,339000,5083000,5086000)
SMODcrop17 <- crop(SMOD_17M, cropbox2)
SMODcrop18 <- crop(SMOD_18M, cropbox2)
SMODcrop19 <- crop(SMOD_19M, cropbox2)
SMODcrop20 <- crop(SMOD_20M, cropbox2)
SMODcrop21 <- crop(SMOD_21M, cropbox2)

NDVIcrop17 <- crop(NDVImax.2017, cropbox2)
NDVIcrop18 <- crop(NDVImax.2018, cropbox2)
NDVIcrop19 <- crop(NDVImax.2019, cropbox2)
NDVIcrop20 <- crop(NDVImax.2020, cropbox2)
NDVIcrop21 <- crop(NDVImax.2021, cropbox2)

4.3.4.2 Visualisation du découpage appliqué

Figure 5: Visualisation de la couche SMOD 2017 après découpage selon une zone délimitée

Figure 5: Visualisation de la couche SMOD 2017 après découpage selon une zone délimitée

4.3.5 Préparation du tableau de données NDVI max

Les données de NDVI max sont récupérées pour chaque année et pour chaque pixel.

NDVI_all <- data.frame(NDVI2017=NDVIcrop17[which(!is.na(NDVIcrop17[]))],NDVI2018=NDVIcrop18[which(!is.na(NDVIcrop18[]))],NDVI2019=NDVIcrop19[which(!is.na(NDVIcrop19[]))],NDVI2020=NDVIcrop20[which(!is.na(NDVIcrop20[]))],NDVI2021=NDVIcrop21[which(!is.na(NDVIcrop21[]))])
head(NDVI_all)
NDVI2017 NDVI2018 NDVI2019 NDVI2020 NDVI2021
0.3653489 0.3707226 0.3768148 0.4216587 0.4651088
0.3420864 0.3607130 0.3040986 0.3711887 0.4289954
0.1543539 0.2082854 0.0792078 0.1866200 0.1487487
0.0955285 0.0887893 0.0729308 0.1797512 0.0875441
0.2148243 0.2155243 0.2773324 0.2985634 0.2387886
0.3455674 0.3172919 0.4380009 0.4517799 0.4957235

4.3.6 Préparation du tableau de données SMOD

De même, les données de SMOD sont récupérées pour chaque année et pour chaque pixel.

SMOD_all <- data.frame(SMOD17M = SMODcrop17[which(!is.na(SMODcrop17[]))], SMOD18M = SMODcrop18[which(!is.na(SMODcrop18[]))], SMOD19M = SMODcrop19[which(!is.na(SMODcrop19[]))], SMOD20M = SMODcrop20[which(!is.na(SMODcrop20[]))], SMOD21M = SMODcrop21[which(!is.na(SMODcrop21[]))])
head(SMOD_all)
SMOD17M SMOD18M SMOD19M SMOD20M SMOD21M
145.000 161.5 140.25 162.4375 116
145.000 165.5 143.25 163.7500 116
145.000 162.5 145.75 165.2500 116
149.125 161.0 149.50 166.0000 116
157.375 161.0 154.50 166.0000 116
157.375 160.5 153.25 163.5625 116

4.4 Calcul des pentes de tendance

La suite du projet consiste à calculer les pentes de tendance sur les 5 années à la fois pour les données NDVI max et SMOD. Le jeu de données est très faible (seulement 5 années) mais l’objectif de ce projet était surtout de définir une méthode d’analyse bien que celle ci soit plus adaptée à une série temporelle plus longue.

4.4.1 Calcul des pentes pour le NDVI max

4.4.1.1 Boucle de calcul des pentes et des pvalue associées

Greening <- data.frame(Pente = rep(NA, nrow(NDVI_all)),
                    Pvalue = rep(NA, nrow(NDVI_all)))

for(n in 1:nrow(NDVI_all)){
  
  ndvi_pixel = data.frame(NDVI = as.numeric(NDVI_all[n,]), Year = c(2017:2021))
  Pente <- summary(mblm(NDVI ~ Year, data=ndvi_pixel))$coefficients[2,1]
  Pvalue <- summary(mblm(NDVI ~ Year, data=ndvi_pixel))$coefficients[2,4]
  
  Greening[n,"Pente"] <- Pente
  Greening[n,"Pvalue"] <- Pvalue
}

4.4.1.2 Génération du raster des pentes NDVI max

NDVI_pente <- NDVIcrop17
NDVI_pente[which(!is.na(NDVI_pente[]))]=Greening[,"Pente"]
Figure 6: Raster des pentes associées à l'évolution du NDVI max entre 2017 et 2021

Figure 6: Raster des pentes associées à l’évolution du NDVI max entre 2017 et 2021

Comme attendu, la plupart des pentes sont proches de 0 traduisant une faible évolution des données NDVI max sur les 5 années considérées.

4.4.2 Calcul des pentes pour le SMOD

4.4.2.1 Boucle de calcul des pentes et des pvalue associées

Melting <- data.frame(Pente = rep(NA, nrow(SMOD_all)),
                    Pvalue = rep(NA, nrow(SMOD_all)))

for(n in 1:nrow(SMOD_all)){
  
  SMOD_pixel = data.frame(SMOD = as.numeric(SMOD_all[n,]), Year = c(2017:2021))
  Pente <- summary(mblm(SMOD ~ Year, data=SMOD_pixel))$coefficients[2,1] 
  Pvalue <- summary(mblm(SMOD ~ Year, data=SMOD_pixel))$coefficients[2,4]
  
  Melting[n,"Pente"] <- Pente
  Melting[n,"Pvalue"] <- Pvalue
}

4.4.2.2 Génération du raster des pentes SMOD

SMOD_pente <- SMODcrop17
SMOD_pente[which(!is.na(SMOD_pente[]))]=Melting[,"Pente"]
Figure 7: Raster des pentes associées à l'évolution du SMOD entre 2017 et 2021

Figure 7: Raster des pentes associées à l’évolution du SMOD entre 2017 et 2021

Les valeurs de pente positives sont plutôt liées aux zones de vallée alors que les valeurs négatives sont plutôt liées aux versants traduisant surtout une précocité de la fonte des neiges sur les versants.

4.5 Etude de la corrélation entre l’évolution du NDVI max et l’évolution du SMOD

Un modèle linéaire est appliqué aux deux variables quantitatives.

## 
## Call:
## lm(formula = df$GreeningPente ~ df$MeltingPente)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.123012 -0.005490 -0.000536  0.005438  0.140579 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      2.442e-03  5.543e-05   44.06   <2e-16 ***
## df$MeltingPente -1.808e-04  4.764e-06  -37.95   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.01646 on 89998 degrees of freedom
## Multiple R-squared:  0.01575,    Adjusted R-squared:  0.01574 
## F-statistic:  1440 on 1 and 89998 DF,  p-value: < 2.2e-16
Figure 8: Relation entre NDVI max et SMOD

Figure 8: Relation entre NDVI max et SMOD

Une très faible relation négative attendue est observée entre l’évolution du SMOD et l’évolution du NDVI max.

Cependant, le modèle a été appliqué à l’ensemble des pixels sans tri préalable des pentes en fonction de la p value associée.

La relation est à nouveau testée seulement pour les pixels qui présentent les pvalue les plus faibles (p value < 0.07).

## 
## Call:
## lm(formula = df2$GreeningPente ~ df2$MeltingPente)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.089626 -0.007452 -0.000791  0.008617  0.082001 
## 
## Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       4.308e-03  2.632e-04   16.37   <2e-16 ***
## df2$MeltingPente -3.759e-04  1.919e-05  -19.58   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.01889 on 5802 degrees of freedom
##   (1689 observations effacées parce que manquantes)
## Multiple R-squared:  0.062,  Adjusted R-squared:  0.06184 
## F-statistic: 383.5 on 1 and 5802 DF,  p-value: < 2.2e-16
Figure 9: Relation entre NDVI max et SMOD pour les p values faibles

Figure 9: Relation entre NDVI max et SMOD pour les p values faibles

Une relation négative est observée entre l’évolution du SMOD et le NDVI max (p value < 0.01). Ainsi, des valeurs de SMOD faibles (dates de fonte précoces) seraient liées à des valeurs de NDVI élevées.
Ce phénomène peut être expliqué par des périodes végétatives plus longues permises par une fonte précoce de la neige qui favoriseraient le verdissement de ces zones.
Cependant, le R2 associé à la relation linéaire est très faible (0.06), probablement en raison du nombre de données élevé issues de zones très diversifiées qui peut augmenter la variabilité observée et du faible nombre d’années étudiées qui ne permet pas d’avoir des évolutions significatives. En effet, peu de variations du NDVI max étaient observées sur les 5 années considérées.
Il faudrait appliquer ce modèle à une série temporelle plus longue ou à une zone d’étude plus homogène (exposition, altitude, type d’habitats) pour vérifier l’existence d’une corrélation entre les deux variables.